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We re-examine the charm contribution to atmospheric lepton fluxes in 
the context of perturbative QCD. We include next-to-leading order correc- 
tions and discuss theoretical uncertainties due to the extrapolations of the 
gluon distributions at small-x. We show that the charm contribution to the 
atmospheric muon flux becomes dominant over the conventional contribution 
CN ' from 7r and K decays at the energies of about 10 5 GeV. We compare our 

qJq ■ fluxes with previous calculations. 
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Neutrino and muon fluxes from cosmic ray interactions with the Earth's atmosphere 
^ ■ have been topics of considerable experimental and theoretical interest M. At energies near 
1 GeV, the 1MB |2| , Kamiokande || and Soudan [[| experiments detect an excess of v e rel- 
ative to in the atmospheric neutrinos. Recent results from SuperKamiokande || appear 
to confirm this observation. At these energies, leptonic decays of charged pions and leptonic 
and semileptonic decays of kaons are responsible for the lepton fluxes, the so-called "con- 
ventional" lepton flux. Currently, it is believed that the conventional flux dominates until 
energies of about 10 3 TeV, when the effects of atmospheric charm production and decay 
become important contributions to the lepton fluxes. The issue of where the charm contri- 
butions dominate is of interest, in part, because this is an energy regime accessible to large 
underground experiments |J. Recent results from Frejus [0], Baksan || and other experi- 
ments show an excess relative to the conventional muon flux in the 10 TeV energy range. 
This may be an indication of a charm contribution at lower energies that expected. One of 
the main goals of the neutrino experiments such as AMANDA |T(| , Antares (ll| , Nestor (12 



X 



and at Lake Baikal |13| are searches for muon neutrinos from extragalactic neutrino sources 



for which atmospheric neutrinos and muons present the main background. 

Lepton fluxes from atmospheric charm have been calculated previously [|14]-17] for specific 



models of charm particle production. Here, we calculate the leptonic flux from charm in 
the context of perturbative QCD. We include next-to-leading order radiative corrections 
and we study the importance of small-x behavior of the parton distribution functions. We 
emphasize the uncertainties inherent in the necessary extrapolation of cross sections and 
energy distributions beyond the experimentally measured regime. We use the comparison 
with low-energy charm production data to constraint some of the theoretical uncertainties, 
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such as the charm quark mass and the factorization and renormalization scale dependence. 
We compare our results to the earlier work on the prompt muons from charm including a 
recent calculation [14fl calculated using PYTHIA Monte Carlo program 



In the next section, we describe the framework for the calculation of the lepton fluxes. 
In Section III, we focus on the charmed quark contribution. In Section IV, we present our 
results for the fluxes and compare with other calculations. We conclude in Section V. 



II. LEPTON FLUX CALCULATION 

Particle fluxes are determined by solving the coupled differential equations that account 
for production, decays and interactions of the particles. The general form of the cascade 
equations describing the propagation of particle j through column depth X is given by 



dX Xj A (dec) 



+ (2.1) 



j 

where Xj is the interaction length, X^^ ~ ^CTjp(X) is the decay length, accounting for time 
dilation factor 7 and expressed in terms of g/cm 2 units. The density of the atmosphere is 
p(X) and 



<t> h (E h ,X)dn k -+j(E;E k 



Je X k {E k ) dE 

In Eq. (2.2), dn/dE refers to either the production distribution l/a k -da k ^j /dE or decay 
distribution l/T k -dT k ^j/dE (where X k — > X k dec ^ in Eq. (2.2)) as a function of the energy E 
of the outgoing particle j. 

It is possible to solve these equations numerically, however, it has been shown JHJ that 
the same results can be obtained with an analytic solution which was derived by noticing 
that the energy dependence of the fluxes approximately factorizes from the X dependence. 
Consequently, one can rewrite 

qfi „ _ ME,X) foo MEk,0) AfcCgj dn k ^(E;E k ) 

b{ ^-^(eTJe dLk ME,0) h(E k ) dE (2J) 

-^ X) -Z kA E). 



X k (E) 

It is often convenient to write Z k j in terms of an integral over xe = E/E k , so 

Z (E)= f 1 dx EME/x E ,0) XkjE) dn k ^(E/x E ) , . 

kA ' Jo x E <j) k (E,0) Xk(E/x E ) dx E 1 ' ] 

In the limits where the flux has a single power law energy behavior, the interaction lengths 
are energy independent and the differential distribution is scaling (energy independent), 
the Z-moment Z k j(E) is independent of energy. In practice, the Z-moments have a weak 
energy dependence because dn/dxE depends on E k , the interaction lengths A are not energy 
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independent, and in general, <pk{E) is not a constant power law in energy over the full energy 
range. The cosmic ray flux can be represented by the following flux of primary nucleons at 
X = 0: 



,{E,X = 0)[cm 



~ 2 s 1 sr 



iGeV- 1 ] = 1.7 (£/GeV)- 2 ' 7 
174 (£/GeV)- 3 



E < Eq 
E>E , 



(2.5) 



where Eq — 5 ■ 10 6 GeV pT| , |2"2"|| . At these energies, we assume isotropy of the flux f23| . 

The detailed solutions to the cascade equations can be found, for example, in Refs. 
|T9| and POfl . Following Ref. [14], we assume that the incident cosmic ray flux can be 
represented by protons. The flux results, in high energy and low energy regimes for lepton 
flavor £ = u^, u e or p due to proton production of hadron j followed by j decay into £ are 



ij,high Z p j {E)Z je {E) Mkj/KJmjcho 

* l-Z pp {E) l-V A i Er J 



f(9) rp (E,0), 



(2.6) 



4 ^ = Z^(E)^) MEtQh (27) 

where an isothermal model for the atmosphere, in which p(h) = poexp(—h/h ) describes 
the density profile as a function of altitude h. The parameters are h = 6.4 km and po 
'2.0'.] x 10~ 3 g/cm 3 p4 |. The quantity rrij is the decaying particle's mass and 



is an effective interaction length, which is weakly dependent on energy. The zenith angle 
dependence of the high energy flux is characterized by f(6) ~ 1/ cos8 for 9 < 60°. At 
higher zenith angles, f(9) is a more complicated function which accounts for the curvature 
of the earth. Details appear in Ref. p0| . The low energy flux is isotropic. When the 
cascade involves charmed hadrons, the low energy behavior dominates and the flux is called 
"prompt" . Critical energies, below which the decay length is less than the vertical depth of 
the atmosphere, range from 3.7 — 9.5 x 10 7 GeV [14]. Interpolation between high and low 
energy fluxes is done via 

0Ji low high 

^ e ~ ajJow i ijMgh ' (2-9) 
j H>i t re 

Eqs. (2.6) and (2.7) show that the bases for the calculation of the prompt lepton fluxes 
are production and decay Z-moments involving charm. The decay moments are discussed in 
Section III.D. The main uncertainties in the calculation of the lepton flux from atmospheric 
charm are the production Z-moments: Z p £, and Z p \ c . The production moments are given 
by 

z =2 f 1 dxE<t>p{E/x E ) 1 da pA ^ c5 {E/x E ) 
PC Jo x E (j)p(E) cr pA (E) dx E 
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The differential cross section is evaluated here using perturbative QCD. The factor of two 
accounts for the multiplicity of charmed (or anticharmed) particles. The charm Z-moments 
can be converted to hadronic moments by 

Z pj (E) = ft Z pc (E) , (2.11) 

where fj is the fraction of charmed particles which emerges as hadron j, where j = 
D°, D + , and A c . We implicitly sum over particles and antiparticles (hence the fac- 
tor of two in Eq. (2.10)). 

The inelastic proton- air cross section <j p a{E) is parameterized by [^5| 

a pA (E) = 280 - 8.71n(£/GeV) + 1.14 ln 2 (£/GeV) mb . (2.12) 

In the high energy limit of the lepton fluxes, in addition to Z pc , we need effective hadronic 
interaction lengths Aj. The proton effective interaction length is therefore 

K ' {E) * iCwa\) (2 ' 13) 

where A = 14.5 is the average atomic number of air nuclei and Nq = 6.022 x 10 23 /g. We use 
Thunman et a/.'s (TIG) energy dependent Z pp , calculated using a PYTHIA Monte Carlo 
|1~8| as a function of energy [[14]]. The charmed hadron j interaction lengths are all taken to 



be equal to the kaon interaction length, approximated by 

A e pp {E) 1 

A,' ^ ; : Z, r . 2.14 

3 N a pA (E)a% P (E)(l-Z KK ) 1 ; 



We use Zkk from Ref. |14]]. The total cross sections are parameterized using the particle 
data book values JMJ based on Regge theory [p7[ . The prompt lepton flux below 10 8 GeV is 
insensitive to the detailed values of Aj because essentially all of the charmed hadrons decay 
before reaching the surface of the earth. Therefore, for most of the energy range considered 
here, the charmed particles are "low energy" and Eq. (2.7) describes the lepton fluxes. 
We now turn to the evaluation of Z pc in perturbative QCD and the other charm inputs. 



III. CHARM CROSS SECTION AND ENERGY DISTRIBUTION 

The charm production cross section and energy distribution are the largest uncertainties 
in the calculation of the prompt lepton fluxes. Since the charm quark mass is of the order 
of 1.3 GeV, the treatment of the charm quark as a heavy quark may be questionable. 
Theoretical uncertainties, due to the possible range of charm quark masses, as well as the 
usual factorization and renormalization scale dependence need to be studied. Theoretical 
predictions based on perturbative QCD calculation fit the available data reasonably well 
in the energy range up to 800 GeV beam energy PBfl . However, atmospheric lepton flux 



calculations require beam energies up to and beyond 10 s GeV. The parton distribution 
functions are needed at very small parton momentum fraction x, outside of the measured 
regime 



29 



In this section we will address these theoretical issues: 
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• the effect of next-to-leading order corrections on the cross section and charmed particle 
energy distribution, 

• charmed quark mass dependence, 

• factorization and renormalization scale dependence, 

• the consequences of the small-x behavior of the parton distribution functions on the 
interaction Z pc moment, and 

• the A dependence of the proton-air charm production cross section. 

From our evaluation of these quantities, a theoretical uncertainty associated with perturba- 
tive charm production will be evaluated. We also describe our inputs to the decay moments 
of charmed hadrons. 



A. Total Cross Section 

The next-to-leading order (NLO) total charm cross section has been calculated by Nason, 
Dawson and Ellis [BO] and by van Neerven and collaborators BTfl. The NLO cross section 



is a factor of between 2 and 2.5 larger than the leading order cross section. Gluon fusion 
dominates the production process. In Fig. 1, we show the importance of the charm quark 
mass in the NLO cross section. We compare the NLO a(pN — » ccX) as a function of the 
beam energy E obtained with the renormalization scale /i equal to the factorization scale M 
equal to the charm quark mass m c with m c = 1.3 GeV and m c = 1.5 GeV. The cross sections 
are evaluated using the CTEQ3 parton distribution functions p2| . The corresponding value 



of Af 5 is 239 MeV. Fixed target data from a summary by Frixione et al. f33[ are also 
plotted. We note that the fixed target data seem to prefer m c = 1.3 GeV. In all of the 
subsequent figures, we set m c = 1.3 GeV. The CTEQ3 parton distribution functions will 
be our canonical set, in part because they incorporate global fits to HERA data, and while 
their validity is not claimed for parton fraction x below x mi „ = 1CT 5 and Qq = 1.6 GeV, the 
program nevertheless provides smooth parton distribution functions below these values. 

In Fig. 2 we show dependence of the total cross section on the scale and parton distri- 
bution. We plot the NLO cross section for different values of /i and M: using the CTEQ3 
structure functions, we set /i = M = m c (dot-dashed) and \i = m c , M = 2m c (solid) with 
m c = 1.3 GeV. The dashed line is the cross section obtained with the MRSD- parton dis- 



tribution functions |34j] and scales /i = m c , M = 2m c with m c = 1.3 GeV. Also plotted are 
the data as in Fig. 1. 

The MRSD- distribution functions have a small-x behavior that is suggested by the 
BFKL approach P5 |. In the small-x limit, the parameterization of the gluon (and sea 



quark) distribution functions at reference scale Qo is 

xg(x,Q ) ~ aT A . (3.1) 



The D- distributions have A = 0.5. Typically, global fits such as the MRSA [§§, MRSG |37| 
and CTEQ3 distributions have A ~ 0.3. By using the D- distributions, we are effectively 
setting an upper limit on the perturbative charm cross section, given our choices of m c , /i 
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and M. We note that, generally, parton distribution functions begin evolution at Qo larger 
than 1.3 GeV. Consequently, our default factorization scale is M = 2m c so that we can use 
more than the CTEQ3 parameterizations. 

Fig. 2 indicates that a low energies, the total cross section has weak dependence on the 
choice of the scale and the parton distribution function. At high energies, E > 10 6 GeV, 
there is a factor of 1.7-2.1 increase from M = m c to M = 2m c . The D- cross section is a 
factor of 1.3 larger than the CTEQ3 cross section at E — 10 6 GeV, both with M = 2m c . 
The D- cross section increases more rapidly because of the steeper small- x behavior of the 
parton distribution function and is enhanced by a factor of 2.6 at 10 8 GeV. This gives an 
overall uncertainty of factor of 5.5 at the highest energy of 10 8 GeV. The MRSA and MRSG 
cross sections for M = 2m c lie between the upper and lower curves in Fig. 2. 

The total charm cross section in p-Air collisions, a p A^ C c{E) , can be written as 

OpA^cc = A J (TpN-^cc (3-2) 

We have evaluated the A dependence for charm pair production using a Glauber-Gribov 
model of nuclear shadowing ||38|| . We find that over an energy range of 10 2 — 10 6 GeV, 



7 = 1.0 — 0.8. Since A = 14.5, the shadowing effect is small, so we set 7 = 1. This is 



consistent with recent measurements at E = 800 GeV ||39|| . Low energy measurements at 



larger xe HH indicate smaller 7 values (7 ~ 0.75), which would reduce our flux predictions 
by an overall factor of 0.5. 

We have used a comparison between data and theory for the total cross section to show 
that m c = 1.3 GeV is a reasonable choice, and to estimate the range of cross sections, related 
to the approximate uncertainty in the flux. To evaluate Z pc , we need the energy distribution 
of the charmed particle. In the next section, we discuss the energy distribution of charm 
quarks in NLO QCD. 



B. Charm Energy Distribution 

NLO single differential distributions in charm quark production have been evaluated Na- 



son, Dawson and Ellis |4]J and incorporated into a computer program, which also calculates 
double differential distributions, by Mangano, Nason and Ridolfi [fl2]| . The program is time 
consuming, so we have incorporated NLO corrections to da/dxE by rescaling the leading 
order distribution. The xe distributions at next-to-leading order are well fit by a i^-factor 
rescaling which is a function of xe, where K is defined by 

= da(NLO)/dx E 
X da("LO")/dx E 1 1 

where "LO" means taking the leading order matrix element squared, but using the two-loop 
a s (n 2 ) and the NLO parton distribution functions. K defined this way shows the effects of 
the NLO matrix element corrections. 

Using the NLO computer program with the CTEQ3 parton distribution functions, we 
show our results for K(E,x E ) for E = 10 3 and 2K(E,x E ) for E = 10 6 GeV in Fig. 3. We 
find that K evaluated using the D- and MRSA distributions agree well with Fig. 3. The 
error bars indicate the numerical errors associated with the Monte Carlo integration in the 
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NLO program. At higher energies, the errors become larger for comparable xe because the 
cross section is dominated by small xe- K can be parameterized as 



K(E, x E ) = 1-36 + 0.42 ln(ln(£/GeV)) (3.4) 



-(3.40 + 18.7(£/GeV)-°- 43 - 0.079 ln(£/GeV) 



™1.5 

X TP 



for /i = m c and M = 2m c . The parameterization is shown by the curves in the figure. 

Using the xe and energy dependent if-factor, we plot the charm quark xe distribution 
for E = 10 3 GeV, 10 6 GeV and 10 9 GeV in Fig. 4. The distributions fall rapidly with x E - 
The convolution of the differential distribution with the ratio of proton fluxes and interaction 
lengths, integrated over xe at fixed outgoing charm quark energy, is what is required for the 
Z- moment. 

Fig. 5 shows the differential moment, dZ pc /dxE- The solid lines are for CTEQ3 dis- 
tributions at outgoing charm energies E = 10 4 GeV, 10 6 GeV and 10 8 GeV, in increasing 
magnitudes. The dashed lines represent the same quantities for the D- calculation. The D- 
distributions have approximately the same shape as the CTEQ3 distributions, but there is 
a more rapid growth in overall normalization with energy. 

In the context of perturbative hard scattering production of charm pairs, the average xe 
value in the evaluation of Z pc is 0.15-0.2. More than 80% of the cross section comes from 
charm transverse momenta below a value of 2m c . In the low transverse momentum limit, 
xe — xe- Fixed target experiments measure da/dxF- The measured charmed meson xf 
distributions are consistent with the perturbative NLO QCD calculations for charm quark 
production, without any fragmentation corrections that would soften the xe distributions 
f43fl . Fragmentation calculations are applicable at large transverse momentum. For the 
calculation of Z pc , we are in the low transverse momentum regime, so we do not need 
fragmentation. 



C. Hadron Fractions 

We account for the transformation of charmed quarks into hadrons by an energy inde- 
pendent hadronic fractions. The hadronic fractions convert Z pc into the interaction moments 
for the charmed mesons and the A c via Eq. (2.11). The hadron fractions can be obtained 
by the observation that p3| , ^4] 



a(D s ) ~ 0.2 a(D° + D + ) (3.5) 
cr(A c ) ~ 0.3 a(D° + D + ) . (3.6) 

The fractions of charmed quarks that appear as D s and A c are 

f D , = 0.13 , (3.7) 
f Ac = 0.20 . (3.8) 

To get the D + and D° fractions, one needs a ratio of a(D + )/a(D°). Using arguments based 
on isospin invariance and counting of states in the production of D and D*, together with 
branching fractions for D* — > D, Frixione et al. suggest that <j(D + )/<j(D°) ~ 0.32. 



With this assumption for the ratio of the cross sections, 
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/do = 0.51 , 
f D+ = 0.16 . 



(3.9) 
(3.10) 



There is some uncertainty in the values of fj. Experimental measurements of 
a(D + )/a(D°) in pN and pp fixed target experiments tend to lead to a somewhat higher 
ratio of cross sections. For example in pp collisions with a beam energy of E\, = 400 GeV, 
the LEBC-EHS Collaboration @| measures a(D+)/a(D°) = 0.7 ± 0.1, while in pN colli- 
sions at Efj = 250 GeV, the ratio is measured by E769 K| to be 0.57 ± 0.22. By taking 
a(D + )/a(D°) = 0.6, the resulting change in the predicted flux is only ~ 15%. 

Integrated Z-moment Z p e>o scaled by 10 3 versus charmed particle energy is shown in 
Fig. 6 for the D- and CTEQ3 distributions with /x = m c and M = 2m c . Also shown is the 
CTEQ3 calculation with fi = M = m c . The other Z-moments for charm production are 
simple rescalings of Fig. 6. While the curves are similar up to 1 TeV, by E = 10 6 GeV, 
there is a factor of ~ 5 between the upper and lower curves. At E = 10 8 GeV, the upper 
and lower curves differ by more than a factor of 10, larger than the ratio of cross sections at 
the same energy. This is accounted for by the fact that the Z-moments at energy E involve 
integrals of the cross section at a higher energy. In addition, since xe ~ 0.15 — 0.2 and the 
cross section is dominated by parton invariant masses near m c , small parton x values are 
emphasized. Since the prompt flux is proportional to Z pc , this enhancement is reflected in 
the flux as well. 



D. Decay Moments 

The last elements of the calculation of the lepton fluxes from charm are the decay mo- 
ments Zki(E) for k = D + , D°, and A c . The decay moments can be written in the same 
form as Eq. (2.4) with X k now representing the decay length. The decay distribution can 
be represented by 



dn k ^e(E; E k ) _ 
dE 

so the decay moments, in terms of an integral over xe = E/E k , are 




(3.11) 



Z,,iE) = j dx E M f/* E) F k ^(x E ) . (3.12) 



The function F is given in Ref. p0|J17| , in the approximation that the leptonic decays of 
charmed mesons are approximated by three-body decays. In the ultrarelativistic limit, F is 
nearly equal for £ = v e , and /i, so the decay moments for the three leptons are essentially 
equal. Consequently, we take the fluxes for the three leptons to be equal. 



Following Bugaev et al. in Ref. [0, the effective hadronic invariant mass mx of the decay 
of the D + is taken to be 630 MeV and for the D° decays, 670 MeV. We take m x = 840 
MeV for D s decays. The A c will ultimately contribute little over most of the energy range 
considered. For A c decays, we use the same three-body formula with an effective hadronic 
mass of 1.3 GeV. The branching ratios are: 
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B(D + 


= 17% 


B(D° -> 


£) = 6.8% 




= 5.2% 




t) = 4.5% 



(3.13) 



For the energies considered here, the charmed particle fluxes in the low energy limit 
dominate. Assuming that Z pp and A p are nearly energy independent, this means that 

<Pk ~ Z pk (E) E(f) p (E) . (3.14) 

The proton flux falls like E~ 2 — E~ 3 . The charm production Z-moments increase with 
energy, as seen in Fig. 6. When we put in the low energy D + meson flux and evaluate the 
Zp+i moment, we get the results shown in Fig. 7. All of the other low energy decay moments 
can be obtained by branching fraction rescaling. For the high energy moments, we take 
Z pk ~ E 0A2 for the D- distributions and Z pk ~ E - 23 for CTEQ3, with fc ~ Z ph {E)(j) p {E). 



IV. PROMPT LEPTON FLUX 

In Fig. 8 we show our results for the prompt atmospheric flux scaled by E 3 for two 
parton distributions and factorization scale choices. The highest flux at E — 10 8 GeV is 
with the D- distribution and M — 2fi — 2m c (dashed). The CTEQ3 distributions with the 
same choice of scale are represented by the solid line, while the dot-dashed line shows the 
result when M — \l = m c . For reference, we show the vertical conventional and prompt flux 



calculated and parameterized by TIG in Ref. [|14|| . The fluxes directly reflect the interaction 



Z-moments of Fig. 6. We emphasize that the prompt flux is isotropic except at the highest 
energies, while the conventional flux is not. 

We have also estimated the flux due to pion-air interactions creating charm pairs. The 
effect is to increase the prompt flux by ~ 30% at 10 2 GeV and by ~ 15% at 10 6 GeV. This 
is a small effect, so we neglect pion contributions to charm production. 

The prompt lepton flux evaluated using perturbative QCD can be parameterized as 

log lo (£ 3 £ (£)/(GeV7cm 2 ssr)) = -A + Bx + C x 2 - D x 3 (4.1) 

where x = log 1Q (E / GeV) . In Table I, we collect the constants for the D- and CTEQ3 fluxes 
exhibited in Fig. 8. 



The TIG flux relies on PYTHIA calculations with m r = 1.35 GeV and the MRSG 37 



parton distribution functions. For E < 1 TeV, the TIG flux is larger than ours because 



of fragmentation effects. Their calculation uses the Lund hadronization model |F7| which 
can give the charmed hadron a larger energy than the charmed quark, an effect which is 
larger for smaller center of mass energies. The TIG prompt flux calculation is lower than 
our calculation for energies E > 1 TeV. At parton x < 10 -4 , the parton distributions in the 
TIG calculation are flattened. For example, the gluon distribution is 

xg(x,Q) ~x-° m (4.2) 

below x = 10 -4 . We note that HERA data at 10~ 5 < x and Qq ~ 2 GeV show no indication 



of flattening 29 . NLO effects in the TIG calculation are accounted for by an energy and 
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xe independent factor of 2. Overall, the net effect is that the TIG Z pc is nearly energy 
independent. 

The prompt muon flux by Volkova et al. [|16|] is larger than our calculated flux. This 
comes in part because of their assumption that da/dxE ~ (1 — xe) 5 /%e, independent of 
center of mass energy and a cross section larger than the perturbative one below E ~ 10 5 
GeV. This xe dependence is harder than the perturbative xe dependence shown in Fig. 4. 

Bugaev et al. (|T7| have presented calculations of the prompt muon flux using two phe- 
nomenological, nonperturbative approaches. One is based on the Recombination Quark- 
Parton Model (RQPM) and the other on the Quark-Gluon String Model (QGSM). The 
QGSM prompt flux is relatively small compared to the RQPM flux, which already affects 
the total atmospheric muon flux at energies of a few tens of TeV. Relative to the RQPM 
calculation and the Volkova et al. results, our D- prompt flux is lower. 

Several experiments show an excess in muon flux above ~ 10 TeV [[^H. Following Rhode 
in Ref. |7]], we plot in Fig. 9 the quantity E 3S5 (p^(E), where M represents the sum of the 
prompt and vertical conventional flux. Also shown are the data from Ref. 0. The energy 



scale factor mostly accounts for the rapidly falling conventional flux ||48|| . When we add the 
prompt fluxes of Fig. 8 to the TIG vertical conventional flux, one sees an enhancement at 
muon energies above 10 5 GeV, at a higher energy than the experimental excess shown by 
data points. 



In Ref. |[49|| , we have shown that it is possible to enhance the prompt flux sufficiently to 
account for some of the observed muon excess at a few TeV. This is accomplished by extrap- 
olating the charm cross section at 1 TeV with a faster growth in energy than predicted by 
perturbative QCD. The xe dependence was taken as da/dxE ~ (1 — xe) a - The inputs are 
consistent with fixed target data below 1 TeV beam energies. We found that the predicted 
prompt flux made significant contributions in the region of the observed excess of muons, 
but it does not fully describe the Frejus data 0. These inputs are not consistent with per- 
turbative QCD. The experimental excess of muons cannot be accounted for by perturbative 
QCD production of charm. 

Another implication of the prompt fluxes calculated here is the ratio of the muon neutrino 
to electron neutrino flux. We define 

Ro = ^- (4.3) 

<Pv e 

Using the TIG parameterization of the vertical 8 = 0° neutrino fluxes, adding conventional 
and prompt contributions, R$ as a function of energy is shown in Fig. 10. At a zenith angle 
of 60°, the conventional flux is a factor of two larger. The ratio Rqo is also shown in Fig. 
10. The quantity R is an interesting diagnostic for the onset of prompt neutrino dominance. 
Unfortunately, at these energies, R is difficult to measure. 



V. CONCLUSIONS 

We find that the perturbative charm contributions to lepton fluxes are significantly larger 
than the recent TIG calculation. The prompt muon flux becomes larger than the conven- 
tional muon flux from pion and kaon decays at energies above ~ 10 5 GeV. We set values of 
the charmed quark mass, renormalization scale and factorization scale by fitting the charm 
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production cross section to low-energy data, then we extrapolate to higher energies. We 
find that the NLO corrections give a correction of more than a factor of two which is weakly 
energy and xe dependent. Nuclear shadowing corrections are small for all energies, due to 
the air nucleus being relatively light. The main uncertainty in the perturbative calculation 
of the prompt flux, given fixed charm mass, factorization scale and renormalization scale, 
is the small- x behavior of the parton distribution functions. Different choices of scales and 
distribution functions, extrapolated to low x with the same power law dependence as for 
x > 10~ 5 , yield as much as a factor of ~ 10 discrepancy in the prompt flux at E — 10 8 GeV. 

If the parton distributions are flattened below some critical value x < x c according to Eq. 
(4.2), our curves in Fig. 8 would overestimate the lepton fluxes. For example, if x c = 10~ 5 , 
the calculated prompt flux at E — 10 6 GeV, using CTEQ3 with M = 2/i = 2m c , is ~ 80% 
of the value shown by the solid curve in Fig. 8, reducing to ~ 40% of the value at 10 s GeV. 
For x c = 10~ 6 , the flux is ~ 70% of our calculated value at 10 8 GeV. However, such an 
abrupt turnover in the power law behavior of the small- a; parton distributions is unlikely. 

We conclude that the prompt muon flux calculated in the context of perturbative QCD 
cannot explain the observed excess of muons in the TeV region |7|-§], independent of the 
theoretical uncertainties associated with small parton x. However, prompt fluxes calculated 
using non-perturbative models of charm production such as discussed in Refs. |TT| , f4"D| ] could 
provide a muon excess in that energy range. Measurements of the atmospheric flux in 
the 100 TeV range would help pin down the charm cross section at energies above those 
currently accessible using accelerators and would provide valuable information about the 
small-x behavior of the gluon distribution function. 

Even though the prompt contributions to the lepton fluxes change the energy behavior of 
the differential fluxes by a factor of E, the atmospheric neutrino fluxes do not compete with 
neutrino fluxes from extragalactic sources above 10 TeV |)TJ. Possible oscillations of muon 
neutrinos as indicated by the Super-Kamiokande experiment || do not affect our results 
due to the extremely small oscillation probability for the energies of interest. 
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TABLES 



TABLE I. 


Parameters for the prompt muon 


plus 


antimuon 


n 

nux appearing 


in big. 8: 


log 10 (£ 3 <^/(GeV 2 /cm 2 s sr)) = -A + Bx + Cx 2 - 


D x 3 , where x = 


log 10 (E/GeV). 




PDF 


Scales A 




B 


C 


D 


CTEQ3 


M = n = m c 5.37 




0.0191 


0.156 


0.0153 


CTEQ3 


M = 2/j, = 2m c 5.79 




0.345 


0.105 


0.0127 


D- 


M = 2/i = 2m c 5.91 




0.290 


0.143 


0.0147 
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FIGURES 
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FIG. 1. The NLO cc production cross section in pN collisions versus beam energy for m c = 1.3 
and 1.5 GeV. The CTEQ3 parton distribution functions are used with M = li = m c . The data are 
taken from the summary in Ref. [28] . 
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FIG. 2. A plot of NLO a^ N versus beam energy for m c = 1.3 GeV using the CTEQ3 (solid) and 
D- (dashed) parton distribution functions with M = 2m c and li = m c . Also show is the CTEQ3 
NLO prediction with M = li = m c (dot-dashed). The data are the same that appear in Fig. 1. 
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CTEQ3, m c =1.3 GeV 
solid: K(10 3 GeV, x E ) 
dashed: 2K(10 6 GeV, x„) 
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FIG. 3. The function K(E,xe) defined in Eq. (3.3) versus xe for E = 10 3 GeV and 10 6 
GeV. The points come from the evaluation of K using the results of Ref. [41,42] with error bars 
indicating numerical errors in the integration, and the curves are our fit to the ratio parameterized 
inEq. (3.4). 




FIG. 4. For E = 10 3 (solid), 10 6 (dashed) and 10 9 GeV (dot-dashed), da/dx El including the 
factor of K(E,xe). The scales used are fi = m c and M = 2m c , for m c =1.3 GeV. 
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FIG. 5. For energies E = 10 4 GeV, 10 6 GeV and 10 8 GeV, dZ pc (E)/dx E versus x E for CTEQ3 
(solid) and D- (dashed) parton distribution functions, where fi = m c and M = 2m c . 
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FIG. 6. Z pD o x 10 3 versus E for CTEQ3 (solid) and D- (dashed) parton distribution functions 
with ix = m c and M = 2m c . Also shown is Z pD o x 10 3 for CTEQ3 with /x = M = m c (dot-dashed). 
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FIG. 7. The decay moment for low energy D + — ► t versus energy for D- (dashed) and CTEQ3 
(solid) with M = 2[i = 2m c and for CTEQ3 with M = \i = m c (dot-dashed). The dotted line 
indicates the decay moment if Z pc is taken independent of energy. 
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FIG. 8. The prompt atmospheric muon flux scaled by versus muon energy for CTEQ3 
(solid) and D- (dashed) with M = 2fi = 2m c . Also shown is the scaled muon flux using CTEQ3 
with M = n = m c (dot-dashed) and the TIG parameterization of the prompt muon flux and the 
vertical conventional muon flux (dotted). 
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FIG. 9. The prompt plus vertical conventional atmospheric muon flux scaled by .E 3 - 65 versus 
muon energy for CTEQ3 (solid) and D- (dashed) with M = 2/j, = 2m c . Also shown is the scaled 
muon flux using CTEQ3 with M = fx = m c (dot-dashed) and the TIG parameterization of vertical 
conventional muon flux (dotted). The data shown are from Ref. [7]. 



upper 0=60 
lower 9 = 




1000 10 4 10 5 10 8 

E [GeV] 

FIG. 10. R = (fru^/tfrve for zenith angles 0° and 60° versus neutrino energy for CTEQ3 (solid) 
and D- (dashed) with M = 2/x = 2m c , and using CTEQ3 with M = \i = m c (dot-dashed). 
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